package org.platform.hrv.lib.common;

import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.TransformType;

import java.util.Arrays;

public class HRVCalculator {
    public static void main(String[] args) {
        // 输入的RRI序列（单位：毫秒）
        double[] rrIntervals = {
                741,753,644,780,732,711,700,680,976,532,829,713,637,406,670,776,814,694,748,462,559,713,472,997,713,713,693,661,745,684,676,710,662,695,695,707,685,681,697,692,721,684,709,720,702,712,695,679,705,713,651,720,734,668,707,955,494,741,733,705,680,693,420,1171,685,528,796,533,1101,1140,845,582,713,527,1102,748,605,665,756,476,650,572,647,713,875,557,713,736,762,749,795,771,778,711,726,761,758,753,727,780,743,742,736,1143,736,686,714,718,745,748,742,787,794,792,730,768,750,755,719,705,740,754,738,718,747,747,750,737,742,740,784,733,703,738,720,731,680,699,718,735,753,745,740,708,710,744,722,724,722,740,745,761,723,717,741,722,705,674,727,707,733,700,730,729,673,770,737,755,745,758,710,731,766,486,670,673,525,673,688,729,710,729
        };

        // 计算时域指标
        double meanNni = calculateMean(rrIntervals);
        double medianNni = calculateMedian(rrIntervals);
        double rangeNni = calculateRange(rrIntervals);
        double sdnn = calculateSDNN(rrIntervals, meanNni);
        double cvnni = sdnn / meanNni;
        double rmssd = calculateRMSSD(rrIntervals);
        double sdsd = calculateSDSD(rrIntervals);
        double[] nn50Results = calculateNN50(rrIntervals);
        double nn50 = nn50Results[0];
        double pnn50 = nn50Results[1];
        double[] nn20Results = calculateNN20(rrIntervals);
        double nn20 = nn20Results[0];
        double pnn20 = nn20Results[1];
        double meanHr = 60000.0 / meanNni;
        double minHr = 60000.0 / Arrays.stream(rrIntervals).max().getAsDouble();
        double maxHr = 60000.0 / Arrays.stream(rrIntervals).min().getAsDouble();
        double stdHr = calculateStdHr(rrIntervals);
        double cvsd = sdsd / meanNni;

        // 计算频域指标 - 使用RR间期序列
        double[] powerResults = calculatePowerSpectrumRR(rrIntervals);
        double totalPower = powerResults[0];
        double vlf = powerResults[1];
        double lf = powerResults[2];
        double hf = powerResults[3];
        double lfHfRatio = lf / hf;
        double lfnu = (lf / (lf + hf)) * 100.0;
        double hfnu = (hf / (lf + hf)) * 100.0;

        // 打印结果
        System.out.println("时域指标:");
        System.out.printf("mean_nni: %.2f ms%n", meanNni);
        System.out.printf("medianNni: %.2f ms%n", medianNni);
        System.out.printf("range_nni: %.2f ms%n", rangeNni);
        System.out.printf("SDNN: %.2f ms%n", sdnn);
        System.out.printf("cvnni: %.6f%n", cvnni);
        System.out.printf("RMSSD: %.2f ms%n", rmssd);
        System.out.printf("SDSD: %.2f ms%n", sdsd);
        System.out.printf("NN50: %.0f%n", nn50);
        System.out.printf("pNN50: %.2f%%%n", pnn50);
        System.out.printf("nni20: %.0f%n", nn20);
        System.out.printf("pnni20: %.2f%%%n", pnn20);
        System.out.printf("mean_hr: %.2f bpm%n", meanHr);
        System.out.printf("min_hr: %.2f bpm%n", minHr);
        System.out.printf("max_hr: %.2f bpm%n", maxHr);
        System.out.printf("std_hr: %.2f bpm%n", stdHr);
        System.out.printf("cvsd: %.15f%n", cvsd);

        System.out.println("\n频域指标:");
        System.out.printf("total_power: %.15f ms²%n", totalPower);
        System.out.printf("vlf: %.15f ms²%n", vlf);
        System.out.printf("LF: %.15f ms²%n", lf);
        System.out.printf("HF: %.15f ms²%n", hf);
        System.out.printf("lf_hf_ratio: %.15f%n", lfHfRatio);
        System.out.printf("lfnu: %.15f%%%n", lfnu);
        System.out.printf("hfnu: %.15f%%%n", hfnu);
    }

    // 计算均值
    private static double calculateMean(double[] values) {
        return Arrays.stream(values).average().orElse(0.0);
    }

    // 计算中位数
    private static double calculateMedian(double[] values) {
        double[] sorted = values.clone();
        Arrays.sort(sorted);
        int n = sorted.length;
        if (n % 2 == 0) {
            return (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0;
        } else {
            return sorted[n / 2];
        }
    }

    // 计算极差
    private static double calculateRange(double[] values) {
        double max = Arrays.stream(values).max().getAsDouble();
        double min = Arrays.stream(values).min().getAsDouble();
        return max - min;
    }

    // 计算标准差 (SDNN)
    private static double calculateSDNN(double[] values, double mean) {
        double sum = 0;
        for (double val : values) {
            sum += Math.pow(val - mean, 2);
        }
        return Math.sqrt(sum / values.length);
    }

    // 计算相邻差值的均方根 (RMSSD)
    private static double calculateRMSSD(double[] values) {
        double sumSquares = 0;
        for (int i = 1; i < values.length; i++) {
            double diff = values[i] - values[i - 1];
            sumSquares += diff * diff;
        }
        return Math.sqrt(sumSquares / (values.length - 1));
    }

    // 计算相邻差值的标准差 (SDSD)
    private static double calculateSDSD(double[] values) {
        double[] diffs = new double[values.length - 1];
        for (int i = 1; i < values.length; i++) {
            diffs[i - 1] = values[i] - values[i - 1];
        }
        double meanDiff = Arrays.stream(diffs).average().orElse(0.0);
        double sum = 0;
        for (double d : diffs) {
            sum += Math.pow(d - meanDiff, 2);
        }
        return Math.sqrt(sum / diffs.length);
    }

    // 计算NN50和pNN50
    private static double[] calculateNN50(double[] values) {
        int count = 0;
        for (int i = 1; i < values.length; i++) {
            if (Math.abs(values[i] - values[i - 1]) > 50) {
                count++;
            }
        }
        double pnn50 = (count * 100.0) / (values.length - 1);
        return new double[]{count, pnn50};
    }

    // 计算NN20和pNN20
    private static double[] calculateNN20(double[] values) {
        int count = 0;
        for (int i = 1; i < values.length; i++) {
            if (Math.abs(values[i] - values[i - 1]) > 20) {
                count++;
            }
        }
        double pnn20 = (count * 100.0) / (values.length - 1);
        return new double[]{count, pnn20};
    }

    // 计算心率标准差
    private static double calculateStdHr(double[] rrIntervals) {
        double[] hrValues = new double[rrIntervals.length];
        for (int i = 0; i < rrIntervals.length; i++) {
            hrValues[i] = 60000.0 / rrIntervals[i]; // 心率 = 60000 / RRI (ms)
        }
        double meanHr = Arrays.stream(hrValues).average().orElse(0.0);
        double sum = 0;
        for (double hr : hrValues) {
            sum += Math.pow(hr - meanHr, 2);
        }
        return Math.sqrt(sum / hrValues.length);
    }

    // 使用RR间期序列计算功率谱
    private static double[] calculatePowerSpectrumRR(double[] rrIntervals) {
        // 步骤1：计算时间点（累积时间）
        double[] times = new double[rrIntervals.length];
        double totalTime = 0;
        for (int i = 0; i < rrIntervals.length; i++) {
            times[i] = totalTime;
            totalTime += rrIntervals[i];
        }

        // 步骤2：重采样设置（4 Hz，即250ms间隔）
        double samplingInterval = 250; // ms
        int numSamples = (int) (totalTime / samplingInterval) + 1;
        double[] resampledRR = new double[numSamples];
        double[] resampledTime = new double[numSamples];
        for (int i = 0; i < numSamples; i++) {
            resampledTime[i] = i * samplingInterval;
        }

        // 步骤3：线性插值RR间期
        int idx = 0;
        for (int i = 0; i < numSamples; i++) {
            double t = resampledTime[i];
            // 找到当前时间点所在的区间
            while (idx < times.length - 1 && times[idx + 1] < t) {
                idx++;
            }
            if (idx == times.length - 1) {
                resampledRR[i] = rrIntervals[rrIntervals.length - 1];
            } else {
                double t0 = times[idx];
                double t1 = times[idx + 1];
                double rr0 = rrIntervals[idx];
                double rr1 = rrIntervals[idx + 1];
                // 线性插值公式
                resampledRR[i] = rr0 + ((rr1 - rr0) * (t - t0)) / (t1 - t0);
            }
        }

        // 步骤4：去线性趋势
        double trendSlope = 0;
        double trendIntercept = 0;
        double meanResampledTime = Arrays.stream(resampledTime).average().orElse(0);
        double meanResampledRR = Arrays.stream(resampledRR).average().orElse(0);
        double cov = 0, var = 0;
        for (int i = 0; i < resampledTime.length; i++) {
            cov += (resampledTime[i] - meanResampledTime) * (resampledRR[i] - meanResampledRR);
            var += Math.pow(resampledTime[i] - meanResampledTime, 2);
        }
        trendSlope = cov / var;
        trendIntercept = meanResampledRR - trendSlope * meanResampledTime;
        for (int i = 0; i < resampledRR.length; i++) {
            resampledRR[i] -= (trendSlope * resampledTime[i] + trendIntercept);
        }

        // 步骤5：加汉宁窗
        double windowSum = 0;
        for (int i = 0; i < numSamples; i++) {
            double window = 0.5 - 0.5 * Math.cos(2 * Math.PI * i / (numSamples - 1));
            resampledRR[i] *= window;
            windowSum += window * window;
        }

        // 步骤6：使用Apache Commons Math进行FFT
        int n = nextPowerOfTwo(numSamples);
        double[] fftInput = new double[n];
        System.arraycopy(resampledRR, 0, fftInput, 0, numSamples);
        
        FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
        Complex[] fftOutput = fft.transform(fftInput, TransformType.FORWARD);

        // 步骤7：计算功率谱密度
        double[] powerSpectrum = new double[n / 2];
        for (int k = 0; k < n / 2; k++) {
            double real = fftOutput[k].getReal();
            double imag = fftOutput[k].getImaginary();
            double magnitudeSquared = real * real + imag * imag;
            
            // 归一化处理
            powerSpectrum[k] = magnitudeSquared / (windowSum * n);
            
            // 补偿负频率分量（直流分量除外）
            if (k > 0) {
                powerSpectrum[k] *= 2;
            }
        }

        // 步骤8：频段积分
        double fs = 1000.0 / samplingInterval; // 采样频率 (Hz) = 4 Hz
        double df = fs / n; // 频率分辨率
        double vlf = 0, lf = 0, hf = 0, total = 0;
        
        // 积分频率范围 (0.0033 Hz 到 0.4 Hz)
        for (int k = 1; k < n / 2; k++) { // 从1开始，忽略直流分量
            double freq = k * df;
            if (freq < 0.0033) continue; // 忽略超低频
            
            if (freq <= 0.4) {
                double power = powerSpectrum[k] * df; // 转换为绝对功率
                total += power;
                
                if (freq < 0.04) {
                    vlf += power;
                } else if (freq < 0.15) {
                    lf += power;
                } else {
                    hf += power;
                }
            } else {
                break;
            }
        }

        return new double[]{total, vlf, lf, hf};
    }

    // 计算大于等于输入值的最小2的幂次
    private static int nextPowerOfTwo(int n) {
        int power = 1;
        while (power < n) {
            power *= 2;
        }
        return power;
    }
}